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ABSTRACT 

Errors in numerical simulations of gravitating systems can be magnified exponentially 
over short periods of time. Numerical shadowing provides a way of demonstrating 
that the dynamics represented by numerical simulations are representative of true 
dynamics. Using the Sitnikov Problem as an example, it is demonstrated that unstable 
orbits of the 3-body problem can be shadowed for long periods of time. In addition, it 
is shown that the stretching of phase space near escape and capture regions is a cause 
for the failure of the shadowing refinement procedure. 
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1 INTRODUCTION 

The sensitivity which A'^-body integrations exhibit to small 
changes in initial conditions and to numerical errors has 
been an acti v e are a of research since Miller's landmark 
study. iMilleiJ (| 19641 ) demonstrated the exponential diver- 
gence of near-by orbits for systems with A*' ^ 32 and found 
that the separation of nearby orbits increases rapidly when 
close binary interactions occur. He suggested that the diver- 
gence of near-by orbits is too rapid to be solely accounted by 
binary interactions and suggests that there must be a col- 
lective effect to account for the results. However, IStandishI 
([1963) showed that the divergence rate was reduced if the 
potential was replaced with a softened potential and con- 
cluded that the divergence is mainly due to close binary 
interactions. 

The dramatic effects of numerical errors on A'^-body in- 
tegrat i ons w as also demonstrated in an important paper by 
iLecarl l|l968l ). After coordinating a study with 11 different 
integrations of the same 25-body problem for 2.5 crossing 
times, Lecar found that quantities such as half mass radius 
and the moment of inertia can change by as much as 100 
percent. In a study with N = 3, iDeionghe and HutI (|200ll ) 
demonstrated that the amplification of initial errors can in- 
crease by as much as 10^". In addition, they showed that 
the growth of errors during close encounters can be ampli- 
fied by as much as 10**, however some of the growth can be 
recovered after the encounter is over. 

The sensitivity to small changes in initial condi- 
tions and numerical errors is a property associated with 
chaotic systems. A measure of the sensitivity of numer- 
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ical errors can be determined by the Lyapunov expo- 
nent A. Early work suggested that the Lyapunov expo- 
nent is inversely proporti o nal to the c r ossing time tcr 
llKandrup and Smitbl Il99ll: iHeggij Il99ll : iGoodman erall 
Il993r ). However, iGoodman et al.l (1993) suggest a depen- 
dence on A*' of the form A"'^ = tcr/ log A'^ or perhaps 
A~^ = tcr/log(log(A'^)), implying that as A'^ increases the 
rate of separation decreases and the Lyapunov exponent in- 
creases. Tlwj£g{iV}_de2endence was la ter numerically veri- 
fied bv lHemsendorf and MerrittI (Il99ll ). 

Despite the difficulty calculating solutions to A''-body 
integrations, computers still remain a useful tool to study 
self gravitating systems. If numerical errors in numerical so- 
lutions to the A'-body problem cause such drastic changes 
in the actual positions and velocities of particles how can we 
trust the dynamics that these solutions represent? Shadow- 
ing is a way of proving that a true solution to a dynamical 
system follows close to a numerical solution. If true orbits 
can be found close to numerical orbits then the dynamics 
represented by the numerical solutions represents true dy- 
namics. 

This study will discuss the existence of shadow orbits 
for the gravitational 3-body problem. First, definitions and 
concepts related to shadowing of dynamical systems will be 
introduced. Next, a refinement procedure which makes cor- 
rections to numerical orbits to reduce the errors incurred at 
each time step will be presented. The Sitnikov problem will 
then be presented and used as a simple model to discuss 
escape and capture of orbits. An approximate Poincare map 
is then presented to model orbits of the Sitnikov problem 
and will be used in conjunction with the refinement proce- 
dure to discuss the validity of numerical solutions by way of 
shadowing. The failure of the refinement procedure to find 
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shadow orbits will then be discussed and regions of phase- 
space where the procedure fails will be delineated. Finally, it 
will be demonstrated that the shadow times for this problem 
can be modeled as a Poisson process. 



2 SHADOWING 

Consider the autonomous ordinary differential equation 

x = /(x), (1) 

where x G R" and / : R" -> R" is a C^ vector field with 
the associated flow represented by t?*. A sequence of points 
{yk}k^o is said to be a pseudo-orbit if there is an associ- 
ated bounded sequence {hk}kLo of positive time such that, 

ahi. I 



\yk^^-r^(yk)\<&, 



(2) 



for A: = 0, 1, ...,M, where 5 > 0. An example of a pseudo- 
orbit is a numerical solution to U]). To show that a pseudo- 
orbit represents some true dynamics for ([T]), it would be 
enough to show that a true orbit follows close to the pseudo- 
orbit. The pseudo-orbit described above is said to be shad- 
owed by a true orbit if there is a sequence of points {xfcj^lo 
and positive times {ifcj^lQ with i9'''(xfc) = x^+i such that 



|Xfc -yfcl < e, 
and 



(3) 



(4) 

for k — 0,1, .., M and small e > 0. The sequence {xfej^Lg 
is known as a shadow-orbit. The shadow-orbit is a true 
solution to lU. 

The first general contributions made on shado wing for 
dynam ical systems were t he shadowing theorems of lAnosovl 
119671) and Bowcnj (|l972l ). Anosov and Bowen considered 
hyperbolic systems and showed that any pseudo-orbit on 
a hyperbolic invariant set has a shadow-orbit. These the- 
orems were gener alized for pseudo-orbits in the vicinity of 
a hyp erbolic set (|Katdll99l1 : lNadzieialll99ll : [Coomes et all 
1 19951 ). For non- hyperbolic systems or for orbits which are 
far from hyperbolic invariant sets these theorems do not 
apply. Shadowing theorems do exist for pseudo-orbits of 
non-hyperbolic systems a nd usually rely on n u merical ver- 
ification of a theorem JCo omes et al. 1994; Chow ct al. 
1989: Chow and Palmeil Il991. : .Chow and Van Vlecb ,1994 : 
Van VlecmT995l '). 



2.1 Refinement procedure 

Sh adowing TV-body simu l ations was first demon strated 
bv iQuinlan and Tremaini (|l992h and iHaved (|200ll ). Both 
these studies considered the refinement procedure found in 
[Crebogi ct al. (1990) to find numerical shadows for the A'^- 
body problem. The refinement procedure is a noise reduc- 
tion technique which can be used to show the existence of 
shadow-orbits. This procedure will be presented for two di- 
mensional dynamical maps, however the procedure can eas- 
ily be adapted for fiow s and has been extend e d to higher 
dimensional systems bv lQuinlan and Tremaind ([1993) • 

Consider the pseudo-orbit {pk}i'Lo of a map f € R^. 
The goal is to find a new less noisy orbit {pfcj^g close to 
the original orbit. Let e*, represent the one step error where 



efc = Pfc - f(pfc_i). 

The refined orbit is constructed by 

Pfc = Pfc +*fc, 



(5) 



(6) 



where ^k is the correction at time step k. Combine equa- 
tions ((5]) and ((6]) to obtain 



f(Pfc-l) - Gfc - f(pk-l), 



(7) 



where pfe = f(pfc_i). Assuming that the correction, $fc, is 
small, expand f(pfc_i) about Pfe_i in a Taylor series to get. 



f(Pfc- 



f(pfc_i) + ifc_i*fc_i, 



(8) 



where Lk is the linearized map at the fcth time step. Sub- 
stitute (O into (O to obtain 



Lk-i^k-i - Sfc. 



(9) 



It is also assumed that the linearized map has an expanding 
direction, u*,, and a contracting direction, Sfe, at each time 
step k. With this assumption, the objective is to find the 
sequences {^k}kLo s^nd {ek}kLo in the coordinates {ufcj^Lg 
and {efcjfelo by 



*fc = QfeUfc + PkSk 

and 

efe = rikUk -l-CfcSfe- 



(10) 



(11) 



The expanding and contracting directions follow the lin- 
earized maps, 



Ufc + l — ifcUfe, 

and 

Sfc + l ~ ifcSfc. 



(12) 



(13) 



For a random |uoj — 1, equation (|12p gives uj, aligned with 
unstable direction at pk after just a few iterations. Starting 
with a random sm and iterating (|13p backwards gives Sk 
aligned with the stable direction at p^ after a few iterations. 
Substitute dlOj and dTT]) into dH) to get. 



ak+iUk+i + Pk+iSk+i = 

Lk{akUk + /SfeSfc) -I- {rjk+iUk+i + Cfe+iSfe+i) 



(14) 



Substituting (I12p and (|13|l into (|14|l yields recursive rela- 
tionships for {ak}k=o ''■iid {Pk}k=o where, 



Ok+i = \LkUk\ak — rik+i, 
Pk+i ~ \LkSk\l3k — Cfc+i- 



(15) 



Equations (|15|l are made computationally stable by calcu- 
lating the coefficients ak starting with om and iterating 
backwards and the coefficients Pk are calculated by choos- 
ing an initial /3o and iterating forwards. The choice of qm 
and Po are arbitrary and are taken to be um = /Jo = 0. Thus 
the sequence of correction coefficients are given by 



OM = 0, Ok = (ak+i + Vki) /\LkUk\, 
/3o = 0, Pk = \LkSk\l3k - Ck+i, 



(16) 



where the values of rjk and (k can be determined directly 
from (O and (fTTj) . 

Once {pk}k=o h8.s been found, the refinement procedure 
can be iterated. Generally, the number of significant digits 
doubles on each iteration of the process. However, cases have 



Shadowing unstable orbits 3 




h = —\z\ — 



Figure 1. The Sitnikov Problem 



been found where the convergence is much slower or does not 
converge. 

The convergence of the refinement proc edure does not 
in itse lf show the existence of a shadow-orbit. iGrebogi et al.l 
|[1990|) provide a containment procedure in two dimen- 
sions which rigorously proves the existence of a shadow- 
orbit. The containment tech nique was la ter extended to 
three dimensional systems bv iHaved l|200ll 'l. A more practi- 
cal approach_Jo£2iigher dimensional systems was developed 
bv ISauer and York c (1991). They showed that for a given 
pseudo-orbit, if the refinement procedure converges - to ma- 
chine precision - and certain quantities of a theorem remain 

bounded, then the pseu do-orbit has a shadow-orbit. 

It has been found (jQuinlan and Tremainelll992l : iHavej 

I2003I ). that one can tell from the convergence of the refine- 
ment procedure alone whether a given pseudo-orbit can be 
shadowed. So, if iterations of the refinement procedure con- 
verge to a new orbit where the one-step errors are the size 
of machine precision, then it is inferred that a shadow-orbit 
exists for the given pseudo-orbit. The new orbit found by 
the refinement procedure is called a numerical shadow. 



3 THE SITNIKOV PROBLEM 

In this study of shadowing for the 3-body problem, a special 
configuration of the restricted 3-body problem known as the 
Sitnikov problem will be considered. The Sitnikov problem 
is the problem of the motion of a mass- less particle, ms, 
on the axis of symmetry of an equal-mass binary (Figure 
[TJ. Following iMoseiJ (Il973t ). units are chosen such that the 
gravitational constant G = 1 and the total mass M = 1. 
Under these conditions, the equation of motion for mz is 
given by 

(17) 



^i^-H 



where z is the position of ms and r the distance from the 
centre of mass to one of the binary masses. The distance r 
can be approximated to first order in the eccentricity, e, by 



1 



(1 — ecost). 



(18) 



Vr^ + z^ 



(19) 



Taking the plane of motion of the binary (z = 0) as a 
Surface Of Section (SOS), consider a map, : (fo,to) — >■ 
(wi,ti), which takes ma from one crossing of the SOS to 
the next. If ms is on the SOS at time io, </> is a map which 
brings vo = i(to) to time ii > to where vi = z(ii) and 
z{t\) = 0. The map has an open domain Do in which every 
point returns to the SOS. As time enters into the problem 
with period 2-n, Do can be considered in polar co-ordinates 
where the radial variable is v and the angular variable is 
given by t. Alternatively, the domain Do can be considered 
on the surface of a cylinder where the initial position on 
the cylinder is defined by to and Eo- Figure [2] (a) shows the 
domain for in cylindrical coordinates. The colour of each 
point represents the number of periods of the binary before 
escape happens. The green regions represent islands of quasi- 
periodic motion. In Figure |3] an example of a quasi-periodic 
orbit which visits the islands of stability in the vicinity of a 
period 7 orbit is provided. 



3.1 An approximate Poincare map 

lUrminskv and Heggid (|2009l ) demonstrated that the 
Poincare map cj) with (|18|) can be approximated by a 
simplectic map tp : {to,Eo) — >■ (ti,-Ei) where 



Ei/2 = Eo 4- acos(to) -I- 6sin(to) 
i,/i = to + 2C(-.Bl/2)-^''^ 



ci/i 
ii 

El 



= ti/2 + 2C(-Si/2)-3/2^ 

= -Ei/2 — a;cos(ti) -I- 6sin(ti), 



(20) 



and a, b and C are constants. The quantities ti/2 and E1/2 
are approximations of the time and energy values of 7113 at 
a local maximum distance from the SOS. It is clear (Figure 
[2] (a)) that the change in energy of ma from one crossing to 
the next is periodic in time and the trigonometric terms in 
(|20|l can be though of as the lowest order in a Fourier ap- 
proximation to this change. The change in time is obtained 
by approximating the motion of ms as Keplerian. The con- 
stants a and b are proportional to the eccentricity of the 
binary whose values can be shown to be. 



0.149e 
0.5075e 



(21) 



and the constant C — tv/{2\/^). 



3.2 Escape and Capture 



Through interactions with the binary as it crosses the SOS, 
ms can gain sufficient energy such that it leaves the SOS 
and does not return. It can be shown that for some positive 
time t* and positive u — {1 — e)/2, if 



Inn 



1 



z(t*y^ + u 



>o. 



(22) 



and the specific energy of 1713 can be defined by 



then \z{t)\ — >■ as i — >■ 00. Setting 2 = in (|22|l gives 
a lower bound on the velocity of orbits which escape on 
the SOS. The solid black region at the top of Figure [2] (a) 
demonstrates how this condition over estimates the escape 
boundary. All energy and time values in this region do not 
return to the SOS. 
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Figure 2. (a) Domain in (to, Eo)-space for ttie map (f) wliere e = 0.61. Eacti point represents an initial condition and the associated colour 
represents the number of periods of the binary before escape. The green regions towards the bottom of the graph represent quasi-periodic 
orbits which remain bound for all time. The solid black region at the top of the image are initial conditions outside of the domain of <j>. 
The escape criterion used is effective in determining escape but crude in approximating the escape boundary on the SOS. (b) Domain in 
(to, i?o)-space for the map (p where e = 0.61. The colour associated with each initial condition represents the number of periods of the 
binary before escape. In (b), the number of periods of the binary is determined by tjv//27r where M is the number of iterations of the 
map. The green regions represent quasi-periodic orbits which do not escape. 



The map, ifi, provides an accurate way of determining 
escape and capture. From equation (|2Up it is found that the 
mapping ip is defined in a region, 



-Eq < — acos(io) 
for 

to G [0, 2tv] 



6sin(io) := dVo, 



(23) 



(24) 



as time enters into the mapping with period 2n. The curve 
9I>o is the escape boundary. Time and energy values above 
900 are said to have escaped. The domain, Do, can be de- 
fined by (|23|l and (|24|l . Initial conditions in Do are mapped 
into the region, ©i, defined by 



E < -acos(t) + 6sin(t) := dVi, 



(25) 



for t £ [0,27r]. Figure 3] shows how the boundaries 92?o and 
dDi intersect. Orbits are mapped from the region under the 
curve dDo to the region under the curve dDi. The region 
Bo = T^oXDi represents energy and time values for which or- 
bits are captured. In the context of the differential equation, 
these are orbits which come from infinity and get captured 



by the binary. Similarly, the initial conditions in the region 
Bi — T)\\Do are energy and time values for which tp is un- 
defined. Again, in the context of the differential equation, 
the region B\ represents orbits which escape from the sys- 
tem. Finally, note that initial conditions for the differential 
equation are such that i > and 2 = on the SOS. So from 
(|19|l . the initial energy can be bounded from below by, 



£>-l/|r(io)|. 



(26) 



forte e [0,27r]. 

Initial conditions are chosen in Dq with (|26|) for e = 0.61 
and plotted in Figure [2] (b) . The colour of each point rep- 
resents the number of periods of the binary determined by 
Im /iir where M is the number of iterations of the map p. 
The green regions represent stable motion whose orbits re- 
main bounded. As energy increases orbits become unstable 
and escape from the system. Notice the similarities between 
Figure [2] (a) and [2] (b) . Both domains have islands repre- 
senting stable orbits as well as large regions representing 
unstable orbits. In Figure [5] an example of a quasi-periodic 
orbit near a period 7 orbit is provided. In addition to the 
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Figure 3. An example of a quasi-periodie orbit near a period 7 
orbit for equation II17I I on the SOS 2 = 0. Initial conditions are 
z(0) = 1.3, 2(0) = 0.0, e = 0.61 and the phase of the binary is .45 
radians from pericentre. 
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Figure 5. An example of a quasi-periodic orbit near a period 
7 orbit for the map ip. Initial conditions are to = 6.01822 and 
Eo = -2.5297 for e = 0.61. 




the refinement procedure can be directly applied as shown 
in section [2?T] 



Figure 4. The curve BTIq represents a lower bound of energy and 
time values for which i/p is undefined. Similarly, the curve dTDx 
represents a lower bound of energy and time values for which the 
inverse map ip"^ is undefined. The shaded region is the domain 
Do for the map 1/3. The two regions labeled Bo and B\ bounded 
by the curves &Do and dTDx axe the capture and escape regions 
respectively. 



similarities between Figure [2] (a) and (^b). IUrminskvl (|2009l ) 
demonstrates t hat th e map ip satisfies Lemmas similar to 
Lemmas 1-5 in iMoserl (|l973l ) (pages 87-91) and that tp, like 
4>, possesses a hyperbolic invariant set on which ip is topo- 
logically equivalent to the shift map. 



4.1 Long lived orbits 

Usin g the co ntainm ent and refinement procedure, 
Gre bogi et al.l (|l990f l successfully demonstrated the 
existence of shadows for pseudo-orbits of length 10 or 
more. To test the algorithm the refinement procedure is 
applied to long lived orbits of the map ip. As seen in Figure 
[2](b), there are regions of stable motion where orbits remain 
bounded forever. The refinement procedure is applied to 
these orbits and it is found that most can be shadowed for 
many iterations. Some of these are sh own in Figure (5] 

As shown bv lPvorak et al.l (|l99a l for the Sitnikov prob- 
lem, the map ip has 'sticky' regions where orbits can be 
trapped for long periods of time before escape. In Figure [S] 
an example of a sticky orbit trapped in the vicinity of is- 
lands of stable quasi-periodic orbits is shown. The inset plot 
in Figure [5] is a magnification of the orbit near one of the 
islands. By sampling the phase space around the islands of 
stable motion, one can find many sticky orbits which survive 
for long periods of time. In Figure [7] the shadow distance is 
plotted against the number of iterations of the map for sev- 
eral sticky orbits where e = 0.61. It is shown that as the 
number of iterations increases, the distance of the numer- 
ical shadow from the pseudo-orbit increases proportionally 
to the number of iterations. 



4 RESULTS 

The map ^p provides a simple way of studying shadowing for 
orbits like those of the Sitnikov problem. The approximate 
map is used to avoid integrating between successive cross- 
ings of the SOS thus obtaining a tremendous speed up in 
calculations. In addition, the one step error can more eas- 
ily be controlled. At each time step uniformly distributed 
noise |(5fe| ^ 5 is added to generate the pseudo-orbit. The 
refinement procedure is then used to reduce the noise level 
to machine precision. Since ^ is a 2-dimensional mapping. 



4.2 Shadowing capture orbits 

Consider uniformly distributed initial values in Bo (Figure 
m for e = 0.25. Initial values are iterated forward for a max- 
imum of 100000 iterations up to the penultimate iteration 
before escaping. For each orbit, the number of iterations, 
M, the orbit was 'shadow-able' for as well as the shadow- 
distance are recorded. The orbits are binned into bins of 
length one iteration and averaged over the bin. The results 
are plotted in Figure [8] where the dots represent the average 
shadow-distance at each iteration of the map. Note that as 
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Figure 6. An example of a 'sticky' orbit wiiicii remains close to 
islands of stable orbits for the map ip for 500, 0000 iterations. The 
inset box is a magnification of the upper most island of stability. 
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Figure 8. Each point is average shadowing distance for the as- 
sociate shadow length M for 10^ initial conditions in ©o where 
e = 0.25. 
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Figure 7. Number of iterations verses shadow distance e. 



M increases, there is increasing variability on the distribu- 
tion of average shadow distances. The data can be fit with 
the curve 7 x 10®M which is similar to the results in Fig- 
ure [7| where the shadow distance is proportional to the orbit 
length. 



5 WHERE DOES SHADOWING FAIL? 

Numerical shadows have been found using the refinement 
procedure for orbits whose length exceeds 10"" iterations for 
the map ip. However, what happens when numerical shad- 
ows are not found? What causes the refinement procedure to 
fail? First, it should be noted that the failure of the refine- 
ment algorithm to converge to a numerical shadow does not 
imply that there is not a shadow-orbit for a given pseudo- 
orbit. A shadow may still exist but t he refinement procedure 
was n ot ab le to co nverge towards it. lQuinlan and Tremaind 
l|l992l ) and lHaved (i2003) found that shadowing breaks down 
during close encounters between particles. This is due to the 
stretching of the velocity subspace during a close encounter. 
In the Sitnikov problem, ma interacts with the binary on 
the SOS and the distance separating ms with the binary 



masses is bounded from below (and above) on the SOS. In 
contrast to the problems discussed in the above mentioned 
studies arbitrarily close encounters do not occur in the Sit- 
nikov problem. However, escape and capture occur during 
close encounters with the binary as 7723 crosses the SOS. 
Near the escape and capture boundaries slight changes in 
the energy of ms as it crosses the SOS can lead to signifi- 
cant changes in the duration of successive crossings of the 
SOS. The map provides a simple way of sampling the phase 
space on the SOS to find regions where shadowing is more 
likely to fail. 

Figure |9] shows shadowing results of 10® initial condi- 
tions. The colour of each point represents either success, yel- 
low, or failure, black, of the refinement procedure. Note that 
only orbits which survived more than three iterations of the 
map are considered. This is because the choice of uo and sm 
would infiuence the results for short lived orbits as (|12p and 
(|13|l may not have had enough time to align uj, and Sk in 
the proper directions. From Figure |9] it can be seen that the 
refinement procedure tends to fail near the escape boundary 
c^I^o- Note also that the refinement procedure fails near the 
boundaries of regions containing orbits which escape after 
three or less iterations. 

The reason the refinement procedure fails in these re- 
gions is that there is a stretching of subspace as orbits near 
the boundary SOq. At a given iteration k, the distance from 
boundary, ^I^o, is given by, 

d^\Ek + acos{tk)+bsm{tk)\. (27) 

From the Jacobian of (1201) it can be shown that 



\Lui 



d 



-5/2 



(28) 



Thus, as d — > 0, the correction coefficients a and /3 go to 
and 00 respectively making it more difficult for the refine- 
ment procedure to converge. 

Figure [10] shows the density of successfully shadowed 
orbits based on the closest approach to the boundary dVo 
for increasing eccentricity values. For each shown eccen- 
tricity value, we select 100000 uniformly distributed ini- 
tial conditions in the region defined by to £ (tt, 27r) and 
Eo G {dVo + 2b sm{to) , dVo) . These boundaries describe 
a band of initial conditions bounded above by the escape 
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Figure 9. The figure shows one million initial contiitions for the 
map if where e = 0.61. The map ip was applieci to each initial con- 
dition 50,000 times or until the resulting orbit escaped. The colour 
associated with each point represents the successful application 
of the refinement algorithm. Black represents initial conditions 
where the refinement procedure failed to converge. Yellow repre- 
sents the successful application of the refinement procedure. Only 
orbits which were longer than three iterations of f are considered. 
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Figure 10. Probability density of close approaches to the escape 
boundary for shadow-able orbits. 



boundary. This band also encompasses the capture region 
Bo • The drop in the density to the right of each curve occurs 
at the distance between the lower boundary curve and the 
escape boundary. Note that the density drops off as initial 
conditions approach the escape boundary. Data was fitted 
using a variable bandwidth kernel density function. 



5.1 Probability of capture 

It was found above that as orbits approach the escape 
boundary the likelihood of an orbit being shadowed de- 
creases. This has an impact on the shadow-ability of orbits 
in the capture region Bq. The capture region area is directly 
proportional to the eccentricity of the binary. As e — > 0, 
the initial conditions in So become pushed up against the 
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Figure 11. Fraction of capture orbits shadow-able using the re- 
finement procedure for increasing eccentricities of the binary. 
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Figure 12. Probability density of shadow durations for the map 
!/3 where e=0.61. The amplitude of the one step noise was set 
at 10~^. For shadow durations T < 5000 the density can be 
approximated by an exponential distribution. For larger T the 
density is inversely proportional to T. 



boundary 9Do. It is expected then that for small eccentric- 
ities, orbits would be less likely to be shadow-able. 

To test this hypothesis, 10^ uniformly distributed ini- 
tial conditions are selected in So and iterated forwards until 
each orbit escapes. This is performed for a variety of eccen- 
tricity values and the fraction of shadow-able orbits in each 
case is determined. The results are shown in Figure fTTI The 
fraction of shadow-able orbits increases as the eccentricity of 
the binary increases. This is because the area of the capture 
region increases proportionally to e. As the area increases, 
initial conditions can be selected at a much further distance 
from the escape boundary making them more likely shadow- 
able. 



5.2 Failure as a stochastic process 

The failure of the refinement procedure can happen at any 
point along the orbit and not necessarily at a close approach 
to the escape boundary. The shadow duration is defined as 
the number of iterations for which a given orbit can be shad- 
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Figure 13. Probability density for the lifetime, Xlft=o *ft' ^"'^ ^'^^ 
map ip where e=0.61. The amplitude of the noise is 10~®. It was 
found that the distribution best fit an exponential distribution. 



owed. For an orbit {{ti, _Bi)}^o the shadow duration, T, can 
take on positive integer values T < M. 

Consider the initial conditions for e = 0.61 shown in 
Figure [9] For each resulting orbit, it is determined how long 
the orbit is shadow-able. Figure [12] provides some informa- 
tion on the distribution of shadow lengths. Initial condi- 
tions are chosen in Do and iterated forwards in time us- 
ing l|20p . Each orbit is iterated for 50,000 iterations or un- 
til the solution escapes. The solid line in Figure [12] repre- 
sents the density of the numerical experiments. The spike 
at 50,000 iterations is mostly due to quasi-periodic orbits 
which remain bounded for all time. As shown in Figure 
1121 the density can be approximated, for small iterations, 
by an exponential density function given by ^exp(— ^a;) for 
^ = 0.0019. The inset graph is a magnification of the density 
for 1000 < Ai < 50, 000. In this range, the density function 
is better represented by the function .025/A'/. 

The map approximates the time between crossings of 
the SOS by considering the motion of ma to be Keplerian. 
Instead of considering the distribution of the shadow dura- 
tion in terms of the number of iterations of ip we can instead 
consider the distribution of shadow times, Im where M is the 
number of iterations of the orbits for which it was shadow- 
able. The solid line in Figure [13] represents the probabil- 
ity density of shadow time for the numerical experiments. 
Again, the data can best be approximated by the exponen- 
tial density function for ^ = .0005. The results found here 
are in agr eement, for s mall shadow durations, with previous 
results bv JHavesI (120031 ') which showed that shadow durations 
for larger A'^-body systems have an exponential distribution 
and can be thought of as a Poisson process. 



6 CONCLUSIONS 

The above res ults confirm , for short lived orbits, previous 
investigations ( JHavesI 120031) that showed numerical shadow 
durations, M, for gravitating systems follow a Poisson pro- 
cess with a exponential density function. The result found 
in this study suggests for longer lived orbits, the density 
function is better approximated by a function proportional 
to 1/M. This may be because the population of longer lived 



orbits tends to be dominated by stable orbits, however this 
has not been investigated. 

In section [5] areas of phase-space where the refinement 
procedure is more likely to fail are characterized. These areas 
are near escape boundaries where there is sufficient stretch- 
ing of phase-space to cause the refinement procedure to fail 
to converge to a less noisy orbit. Interestingly this seems to 
be due to the growth of the variational equations over one 
time step. This does not rule out the failure of the the re- 
finement procedure by the accumulative effect of the growth 
of the variational equatio n associated with large Lyapunov 
exponents as discussed bv lZhu and HavesI ((2009,). 

In Figure [TT] it is demonstrated that as the volume 
of phase-space representing capture orbits decreases, it be- 
comes increasingly difficult to shadow capture orbits. This 
is a result of the distribution of failures of the refinement 
procedure seen in Figure 1101 As the volume of phase-space 
associated with capture decreases, capture orbits get pushed 
up against the boundary 9Do where the one-step growth of 
the variational equations causes the refinement procedure to 
fail. 

Finally, it was found that the shadow distance for an 
orbit is proportional to the number of iterations of the map 
(Figures [7] and |S]). It was noticed that if in addition to ii and 
El , orbits were required to be shadow-able at the half steps 
ti/2 and -Ei/2, then initially shadow-able orbits continued 
to be shadow-able. When shadowing at the half step was 
required, the shadow distance typically increased by about 
a factor of two. 

The Sitnikov problem discussed in this study provides 
a straight forward way of characterizing a domain of initial 
conditions as well as regions of stable and unstable motion. 
Work in progress considers slight changes to the Sitnikov 
problem i n order to st udy shadowing of unstable orbits. For 
example, ISoulis et al.l ([2007) consider slight perturbations 
to the mass and position (away from the z-axis) of ma and 
delineate regions of stable and unstable motion. It would be 
expected that, like the results found in this study, shadow- 
ing with the refinement procedure breaks down near bound- 
aries of escape for unstable orbits. In fact, the break down 
of the refinement procedure near escape boundaries would 
be expected for general 3-body configurations. As solutions 
approach parabolic escape boundaries, an orbit can undergo 
increasingly long ejections from the left-over binary system. 
Small changes in the energy of an orbit in this region can 
cause significant changes in the time of return for the orbit. 
If the refinement procedure could make changes to the orbits 
so as to conserve the energy of the ejected body it might im- 
prove the success rate of the refinement procedu re. Finally, 
the Sitnikov 4-body problem (jSoulis et al.ll2008h provides a 
starting point for examining the relationship between the 
shadowing distance and the number of bodies. Extra bodies 
can be added in circular orbits about the center of mass. 
[Haves! (|2003l ) demonstrates that as the number of moving 
bodies in a fixed potential increases the shadow durations 
decrease. It would be of interest to determine if a similar 
relationship holds for the Sitnikov A^'-body problem. 

It should be stressed again that the failure of the refine- 
ment procedure does not necessarily mean that a shadow 
does not exist for a given pseudo-orbit. It may very well be 
that shadows do exist for orbits in regions where the refine- 
ment procedure fails. We are encouraged that this may be 



the case. Both the Sitnikov problem and the approximate 
Poincare map possess a hyperboUc in varia nt set, A, n e ar the 
escape boundaries (see iMosen ()l973l ) and lUrminskvl ()2009l ) 
respectively). Despite the fact t hat A is near th e bo undary 
&Do, t he shadowing theorems bv lAnosovl (|l967l ) and lBowenI 
l|l972h guarantee that any pseudo-orbit on A has an asso- 
ciated shadow-orbit. This demonstrates that being in the 
vicinity on the escape boundary does not necessarily rule 
out the existence of shadow-orbits. 
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